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ABSTRACT 


An  extensive  review  of  the  literature  revealed  that  many 
experimental  studies  have  been  conducted  in  heat-transfer  wind 
tunnels  simulating  the  cooling  passages  in  turbine  blades. 
However,  very  few  numerical  studies  have  been  performed. 
Phoenics,  a  computational  fluid  dynamics  computer  program, 
produced  results  for  several  duct  configurations  and  calculat¬ 
ed  the  heat  transfer  characteristics  of  each.  The  configura¬ 
tions  investigated  included  a  straight,  square  duct  and  a 
rectangular  duct  with  turbulators  (a  form  of  turbulence 
promoters)  present.  The  parameters  varied  included  Reynolds 
numbers,  turbulence  intensity  and  grid  geometry.  Results  for 
the  turbulated  duct  indicated  highly  distinct  and  repeatable 
flow  patterns  developed  over  a  wide  range  of  values.  The 
variation  of  the  inlet  turbulence  intensity  had  little  impact 
on  the  kinetic  energy.  The  results  proved  to  be  highly  grid 
dependent,  which  greatly  impacted  the  correlation  between 
experimentally  and  numerically  produced  data,  for  the  same 
configurations,  under  similar  operating  parameters. 
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I.  INTRODUCTION 


A.  BACKGROUND  TO  TURBINE  INTERNAL  COOLING 

All  modern  gas  turbine  engines  share  three  basic 
components;  the  compressor,  the  combustion  chamber  and  the 
turbine.  From  the  compressed  mainstream  flow  out  of  the 
combustion  chamber  the  turbine  extracts  the  power  to  operate 
the  compressor  and  to  deliver  the  gases  to  either  a  power 
turbine  for  shaft  work  or  to  a  nozzle  for  thrust.  During 
this  process  losses  occur  in  both  the  compressor  and  the 
turbine  sections.  In  the  compressor  they  serve  to  increase 
the  power  absorbed  while  in  the  turbine  they  act  to  decrease 
the  power  output.  To  overcome  these  losses,  more  air  and 
therefore  fuel  flow  are  required  for  a  required  power  level. 
The  compressor  loss  increases  the  temperature  of  the  flow  out 
of  the  compressor  and  into  the  combustor  and  decreases  the 
amount  of  fuel  which  can  be  added  for  a  given  turbine  inlet 
temperature.  Based  on  cycle  analysis,  the  higher  the 
allowable  temperature  into  the  turbine  the  higher  the  power 
output  of  the  unit.  The  turbine  inlet  temperature  constrains 
the  design  of  all  gas  turbine  engines. 

Thus  for  a  given  engine,  the  amount  of  fuel  which  the  fuel 
control  unit  (FCU)  can  add  is  limited  by  the  ability  of  the 
materials  to  withstand  thermal  loads.  Failures  can  result 
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from  effects  such  as  oxidation,  corrosion  and  erosion,  as  well 
as  creep  (due  to  a  combination  of  thermal  and  centrifugal 
loads)  and  fatigue,  which  ensues  from  both  thermal  and  cyclic 
loads.  The  combination  of  these  effects  determine  the 
operating  life  of  the  turbine.  Two  solutions  to  this  problem 
have  emerged.  One  is  the  improvement  of  materials  for  the 
construction  of  turbine  components,  the  other  is  the  use  of 
compressor  bleed  air  to  cool  the  turbine  components . 

Turbine  blade  cooling  will  be  discussed  in  the  following 
order:  external  cooling  air  over  turbines,  impingement 
cooling,  film  cooling,  transpirational  cooling  and  with  the 
final  focus  on  internal  cooling  using  turbulators . 
Experiments  in  turbine  cooling  began  in  the  1940 's  with  the 
first  units  going  into  operation  around  1960.  Initial  engines 
bled  1.5  to  2  %  of  the  compressed  air  from  the  compressor  to 
cool  the  turbine  200-300°  C  [Ref  1]  .  Modern  aircraft  engines 
can  use  as  much  as  20%  of  the  compressor  discharge  flow  [Ref 
2]  .  This  cooling  air  externally  bathes  the  turbine  blade  with 
initial  rotors  receiving  their  cooling  air  directly  from  the 
compressor  discharge.  Two  shortcomings  associated  with  this 
type  of  cooling  include  the  turbine's  reduced  ability  to  drive 
the  compressor  because  of  the  lower  temperature  of  the  flow  as 
it  enters  the  turbine  and  the  mixing  of  the  cooling  and  the 
mainstream  flow  that  results  in  increased  aerodynamic  losses. 
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Rather  than  just  providing  external  cooling  air, 
individual  cooling  of  the  airfoil  itself  is  another  attempt  to 
elevate  turbine  inlet  temperatures.  Impingement  cooling  at 
the  leading  edge  and  at  midchord  consists  of  placing  a  sheet 
metal  insert  into  the  hollow  turbine  blade  with  appropriately 
placed  venting  holes.  Coolant  flows  through  the  root  of  the 
blade  into  the  insert  and  out  the  impingement  holes.  This 
convective  type  cooling  directs  air  at  the  hottest  spots  on 
the  blade  and  creates  airflow  inside  the  blade.  Often  the 
most  difficult  part  of  the  design  analysis  involves 
determining  the  exact  location  of  the  stagnation  point  or 
hottest  spot  on  the  blade.  After  the  cooling  process  occurs 
the  blade  vents  the  coolant  either  at  the  tip  or  at  the 
trailing  edge,  which  then  joins  the  main  flow. 

Film  cooling  also  uses  air  from  the  compressor  routed  into 
the  root  of  the  blades.  In  film  cooling  this  air  travels 
through  the  blade  and  departs  through  openings  designed  to 
produce  a  blanket  of  cooling  air  to  protect  the  blade  against 
the  mainstream  gases.  The  analysis  on  the  design  of  these 
holes  focuses  on  number,  location,  size,  shape  and  orientation 
angle.  With  the  average  cooling  effectiveness  (ACE)  defined 
as  ; 


(1) 
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where  Ttg  is  the  mainstream  gas  total  temperature,  Tm  is  the 
average  surface  temperature  of  the  airfoil  and  Tc  is  the 
coolant  supply  total  temperature,  engineers  usually  specify 
film  cooling  for  an  ACE  of  0.6.  Film  cooling  is  very 
effective  when  combined  with  other  forms  of  internal 
convective  cooling,  but  the  price  in  total  pressure  losses  for 
putting  such  a  quantity  of  cooling  air  into  the  hotter  main 
gas  flow  is  considerable. 

Another  form  of  cooling  is  transpirational  or  aspiration 
cooling.  This  technique  uses  ■.  porous  material  to  construct 
the  airfoil  and  runs  cooling  airflow  into  and  through  it. 
This  combines  both  convective  and  film  cooling,  but  with  some 
obvious  problems  associated  with  the  strength  of  the  airfoil 
construction.  Designers  using  a  woven  wire  material  achieved 
some  success,  however  difficulties  in  the  shaping  of  the 
airfoils  may  preclude  this  type  of  construction . 

The  final  type  of  cooling  discussed  is  internal  cooling 
passages  with  coolant  run  through  the  blade  and  out  at  the  tip 
and  trailing  edge  of  the  blade.  Only  convective  cooling  is 
used  and  chosen  when  the  design  calls  for  an  ACE  of  less  than 
0.5.  Internal  designers  choose  convective  cooling  because  of 
a  limited  supply  pressure  (a  higher  ACE  would  require  a  higher 
pressure) ,  and  a  higher  effectiveness  level  using  just 
convective  cooling  could  lead  to  a  higher  temperature  gradient 
and  thus  aggravate  the  very  thermal  stress  problem  which 
cooling  tries  to  avoid. 
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B .  EXPERIMENTAL  INVESTIGATIONS 

Investigations  into  internal  passages  for  convective  heat 
transfer  fall  into  tnree  sub-areas;  straight  passages, 
serpentine  passages  with  one  or  more  180°  bends  and  single  and 
multi-pass  channels  undergoing  rotation.  The  investigations 
often  include  analysis  of  different  flow  properties,  channel 
shapes,  aspect  ratios  and  wall  roughness,  which  encompasses 
the  shape,  location,  number  and  orientation  of  turbulence 
inducing  structures.  The  continual  refinement  of  the 
operating  characteristics  of  turbine  engines  has  made  this 
last  characteristic  increasingly  important  as  many  of  the 
other  parameters  are  fixed. 

As  early  as  the  1960's  research  centered  on  the  effects  of 
protrusions  into  the  flow  from  the  walls  of  the  channel . 
Wilkie  [Ref  3]  was  one  of  the  first  to  publish  results  of 
experimentation  on  the  effects  of  rib  roughened  ducts  on  heat 
transfer.  The  interest  in  cooling  operations  in  nuclear 
reactor  fuel  elements  resulted  in  the  findings  he  collected 
and  published.  Through  varia^iors  in  duct  equivalent  diameter 
(d)  ,  rib  height  (e)  and  rib  pitch  or  spacing  (p)  he  found 
correlations  between  the  Stanton  number,  the  Reynolds  number, 
the  p/e  ratio  and  the  e/d  ratio.  The  Stanton  number  is 
defined  as  h/ (upCp)  where  h  is  the  heat  transfer  coefficient, 
u  is  the  velocity  in  the  free  stream  at  the  outer  edge  of  the 
boundary  layer,  p  is  the  fluid  density  and  C  is  the  specific 
heat  at  constant  pressure.  The  results  showed  Stanton  number 
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varied  less  with  Reynolds  number  for  small  values  of  e/d 
( < .  0 1 2 )  ,  and  behaved  much  as  on  a  smooth  surface  at  larger 
values  of  e/d  (>.012)  .  In  addition  he  found  the  value  of  the 
p/e  ratio  giving  maximum  Stanton  numbers  varied  with  e/d  and 
Reynolds  number  but  tended  to  approximately  the  value  six  at 
large  e/d  irrespective  of  Reynolds  number.  Norris  [Ref  4] 
also  detailed  the  effect  of  these  " turbulators"  in  his 
investigation.  After  testing  various  shapes,  he  concluded 
that  although  the  local  heat-transfer  coefficient  varied  in  a 
very  non-uniform  manner  in  the  stream-wise  direction,  there 
was  a  value  of  the  pitch-to-width  ratio  for  the  square 
turbulator  which  gave  an  optimum  value  for  the  average  heat 
transfer  as  well  as  friction  factor.  Increasing  the  ratio 
above  this  level  served  only  to  decrease  these  values.  Also, 
he  developed  a  relationship  between  the  friction  factor  of  his 
test  apparatus  and  the  friction  factor  of  a  smooth  duct,  and 
found  an  optimum  value  existed  for  this  ratio.  Increases 
above  this  value  did  not  give  increased  heat  transfer. 

Using  the  Nusselt  number  as  an  indicator  of  heat  transfer, 
Han,  et .  al .  [Ref  5]  varied  their  rib  angle  (a)  and  spacing  to 
determine  their  effects.  With  a  rib-to-height  ratio  of  ten, 
they  found  the  Nusselt  number  to  be  twice  as  high  with  the 
turbulators  placed  90°  to  the  flow  as  for  the  case  without 
turbulators.  Introducing  an  a  of  60°  increased  the  Nusselt 
number  another  5%  while  an  a  of  45°  delivered  an  additional 
25%  over  the  90°  orientation.  Lau,  et .  al .  [Ref  6]  built  upon 
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this  concept  when  they  examined  fully  developed  flow  in  a 
square  channel  with  turbulence  promoters  on  two  opposed  walls, 
using  turbulators  configured  with  ex's  of  45°,  60°  and  90° , 
lengths  running  from  wall  to  wall,  and  shorter  discrete 
turbulators.  Both  configurations  were  tested,  crossed  or 
parallel  with  respect  to  each  other.  They  found  parallel 
discrete  ribs  with  a's  of  60°  provided  the  highest  ribbed  wall 
heat  transfer,  while  the  same  style  of  ribs  canted  at  30° 
provided  the  lowest  level  of  pressure  drop.  In  general  they 
recommended  using  parallel,  angled  discrete  ribs  for  internal 
cooling  passages  as  they  provided  higher  ribbed  wall  heat 
transfer,  lower  smooth  wall  heat  transfer  and  lower  channel 
pressure  drop.  The  optimum  configuration  for  these  discrete 
ribs  was  a  continually  staggered  row  pattern  of  three  ribs 
followed  by  a  row  of  two  ribs. 

Taslim  and  Spring  [Ref  7]  centered  their  straight  square 
channel  experiments  not  so  much  on  the  placement  of  the 
turbulators,  but  rather  the  geometry  of  them.  They  determined 
that  an  optimum  turbulator  spacing  existed  for  every  given 
turbulator  blockage  ratio  and  aspect  ratio.  The  turbulator 
blockage  ratio  was  defined  as  the  ratio  of  the  turbulator 
height  to  the  hydraulic  diameter  of  a  non-turbulated  section. 
They  also  found  that  the  Nusselt  number  decreased  as  the 
turbulator  aspect  ratio  decreased.  This  became  even  more 
critical  as  the  spacing  dropped  below  the  optimum  value.  They 
also  developed  an  alternatively  shaped  turbulator  which  not 


7 


only  provided  the  same  heat  transfer  performance  as  a  larger 
turbulator,  but  also  enhanced  the  castability  and,  thus, 
reduced  the  probability  of  problems  in  manufacturing.  To  do 
this  they  developed  a  turbulator  with  a  narrow  tip  width  to 
enhance  flow  while  at  the  same  time  widening  the  base  to  ease 
the  problems  of  castability.  The  resulting  shape  is 
trapezoidal  with  slightly  rounded  corners  at  the  tip. 

Prior  to  Boyle's  [Ref  8]  report,  most  data  for  turbulators 
came  from  experiments  using  long  passages  with  measurements 
done  in  regions  where  results  did  not  change  with  distance. 
Boyle  initiated  his  report  with  an  investigation  of  a 
serpentine  channel  without  any  turbulators.  He  researched  the 
effect  of  bends  only  on  heat  transfer  coefficients.  He  found 
that  when  varying  the  outside  edge  of  the  turn  from  square  to 
a  round  geometry  he  did  not  alter  the  heat  transfer  rate.  He 
established  that  the  rate  varied  directly  proportional  to  the 
Reynolds  number  and  the  highest  rate  occurred  at  the  beginning 
of  the  second  leg  after  the  bend.  The  addition  of  turbulators 
markedly  increased  the  heat  transfer  rate.  He  also  found  an 
optimum  height  which  required  the  turbulator  to  extend  into 
the  boundary  layer  buffer  region. 

Chandra,  et .  al .  [Ref  9]  did  additional  work  on  a  multi¬ 
pass  serpentine  channel  investigating  the  effect  of  the 
turbulator  angles  on  the  flow.  Using  a  ratio  of  average 
Sherwood  numbers  normalized  to  Sherwood  numbers  for  fully 
developed  turbulent  flow,  they  equated  this  value  to  the 


ratios  of  average  Nusselt  numbers  normalized  with  respect  to 
Nusselt  numbers  for  fully  developed  turbulent  flow.  Here 
they  define  the  Sherwood  number  (Sh)  as: 

Sh=!b?h  (2) 

VSC 

where  is  the  local  mass  transfer  coefficient  or  the  local 
mass  transfer  rate  per  unit  area  divided  by  the  difference  of 
the  coolant  density  at  the  wall  and  freestream;  Dh  is  the 
hydraulic  diameter,  v  is  the  kinematic  viscosity  and  Sc  is  the 
Schmidt  number  for  the  coolant .  The  Schmidt  number  is  the 
ratio  of  the  kinematic  viscosity  to  the  mass  diffusivity. 
Using  the  above  equation  they  were  able  to  make  a  heat  to  mass 
transfer  analogy.  For  the  straight  channel  they  found 
periodic  Sherwood  numbers  on  the  ribbed  walls  in  the  direction 
of  flow,  while  a's  of  60°  and  45°  produced  higher  Sherwood 
numbers  on  the  outer  walls  than  on  the  inner  ones. 
Turbulators  set  perpendicular  to  the  flow  produced  little 
spanwise  variation  of  heat  transfer,  and  the  slight  spanwise 
variations  they  did  find  diminished  somewhat  after  the  180° 
turns .  In  general  they  found  that  with  respect  to  the  area  of 
the  turn  the  overall  average  Sherwood  number  ratio  was 
independent  of  the  rib  angle  but  tended  to  decrease  slightly 
as  the  Reynolds  number  increased. 

Using  a  three-pass  rib-roughened  channel  Han  and  Zhang 
[Ref  10]  established  the  local  pressure  drop  distribution. 
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With  these  data,  they  were  able  to  determine  tne  local  coolant 
flow  rate  distribution.  Their  results  indicated  that  the 
pressure  dropped  linearly  for  the  first  leg  to  the  turn.  Just 
prior  to  the  turn,  a  slight  rise  followed  by  a  sharp  drop 
occurred  into,  chrough  and  just  after  the  bend.  Flow 
recirculation  caused  a  recovery  of  pressure  during  the  second 
straight  channel,  with  a  similar  pressure  drop  occurring  at 
the  second  turn.  They  also  found  the  turbulators'  biggest 
impact  to  be  evident  in  the  first  straight  channel  prior  to 
the  turn.  The  180°  turns  tended  to  dilute  their  impact  on 
subsequent  legs. 

Chyu  [Ref  11]  did  similar  testing  with  both  two-  and 
three-pass  channels  employing  180°  turns.  The  emphasis  of 
these  experiments  was  on  surface  heat  transfer  coefficient 
versus  pressure  distribution.  With  the  two-pass  channel,  he 
found  that  while  the  pressure  loss  attributable  to  the  turn 
was  significant,  it  was  independent  of  the  Reynolds  number. 
His  findings  of  flow  separation,  reattachment,  impingement  and 
secondary  flows  caused  by  the  180°  turn  are  consistent  with 
the  findings  of  Han  and  Zhang.  He  attributed  the  heat 
transfer  levels  observed  to  be  due  to  these  secondary  flow 
effects.  Also,  even  though  he  found  the  heat  transfer  rates 
of  both  turns  to  be  approximately  equal,  the  total  heat 
transfer  of  a  three  pass  channel  was  about  15%  greater  than 
for  a  two  pass  channel. 
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With  a  three-pass  square  channel,  with  two  180°  turns,  Han 
and  Zhang  [Ref  12]  investigated  the  impact  of  angling 
turbulators  60°  on  opposite  walls.  They  also  opted  for  mass 
transfer  to  be  the  measured  experiment  parameter  and  converted 
to  Sherwood  numbers  to  correlate  the  heat  transfer. 
Configuring  the  ribs  both  crossed  and  parallel  to  the  opposite 
wall,  they  found  a's  of  60°  produced  greater  mass  transfer 
coefficients  than  a's  of  90°,  while  parallel  ribs  generated 
greater  coefficients  than  crossed  ones.  Starting  with  a 
smooth  channel,  they  established  the  value  of  the  second  leg's 
mass  transfer  to  be  60%  greater  than  the  first  leg.  Since  the 
third  leg  closely  approximated  the  second,  they  surmised  that 
the  180°  turn  caused  the  difference.  Finding,  in  general,  the 
60°  parallel  orientation  provided  the  highest  mass  transfer 
coefficients,  this  study  showed  rib  angle  orientation  plus 
sharp  turns  in  the  coolant  path  can  increase  or  decrease  the 
mass  transfer  coefficient.  In  addition,  by  correlating 
previous  data  with  the  results  of  this  experiment,  they 
validated  the  application  of  straight  ribbed  channel  data  to 
the  first  pass  in  multi-pass  ribbed  cooling  channels. 

The  final  area  of  investigation  associated  with  internal 
convective  cooling  of  turbine  blades  relates  to  rotating 
elements.  Clifford  [Ref  13]  used  a  rotating,  three  pass, 
smooth  channel  with  airflow  directed  through  it.  Measuring 
the  heat  transfer  along  the  leading  and  trailing  walls  with 
respect  to  rotation,  he  found  rotation  served  to  both  aid  and 
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inhibit  heat  transfer  with  the  reductions  approaching  30%  on 
the  suction  side  of  the  final  leg  and  the  enhancement  closer 
to  35%  on  the  pressure  side  of  the  initial  leg.  Averaging 
his  data  showed  the  effects  of  rotation  reduced  the  heat 
transfer  to  90%  of  the  statically  measured  values. 

Morris  and  Harasgama  [Ref  14]  performed  similar  tests 
using  axial  airflow  directed  both  radially  inward  and  outward. 
Their  results  showed  outwardly  directed  airflow  promoted 
greater  cooling  on  the  trailing  (pressure)  surface  vice  the 
leading  (suction)  surface.  Inwardly  directed  airflow, 
however,  increased  heat  transfer  on  the  leading  surface.  This 
heat  transfer  was  rotational  Reynolds  number  dependent  as  it 
increased  slightly  as  the  number  was  increased.  Like 
Clifford,  they  also  discovered  instances  of  reduced  heat 
transfer  when  compared  with  the  stationary  tube  condition. 
They  also  found  rotational  buoyancy  tended  to  suppress  heat 
transfer  with  radially  outward  flow.  With  respect  to 
stationary  experiments  they  found  no  areas  of  impaired  heat 
transfer . 

The  findings  of  Kheshi  and  Scriven  [Ref  15]  and  Harasgama 
and  Morris  [Ref  16]  were  consistent.  They  found  that  as  the 
Rossby  number  (which  is  the  ratio  of  the  mean  through-flow 
velocity  of  fluid  along  the  duct  to  a  peripheral  velocity 
formed  by  the  product  of  the  angular  velocity  and  the 
hydraulic  diameter  (Qd/V) )  increased,  the  axial  velocity  on 
one  side  decreased  with  respect  to  the  other  side  and  a  cross 
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channel  flowfield  developed  which  was  perpendicular  to  the 
main  flow.  The  steady  state  for  this  condition  showed  that  a 
two-vortex  pattern  developed.  Increasing  the  Rossby  number 
further  generated  a  four-vortex  flow  pattern  in  the  steady 
state.  This  four-vortex  pattern  required  a  shear-free  flow  to 
develop  on  both  the  top  and  bottom  walls.  The  experiment 
suggested  the  flow  pattern  was  a  function  of  the  aspect  ratio 
of  the  channel,  the  Rossby  number  and  the  Ekman  number.  The 
Ekman  number  is  defined  as  the  ratio  of  the  Rossby  number  to 
the  Reynolds  number  [Ref  17]  .  Harasgama  and  Morris  confirmed 
the  effect  of  heat  transfer  on  both  the  leading  and  trailing 
edges  with  either  outward  or  inward  flow.  In  addition  they 
put  a  value  of  10-20%  on  the  decrease  of  the  heat  transfer  on 
the  leading  side  during  outward  flow.  They  attributed  the 
increase  in  trailing  side  heat  transfer  to  the  strengthening 
of  secondary  flow  as  the  average  Nusselt  number  increased. 

Hajek  et .  al.  [Ref  18]  investigated  the  impact  of  tne  rate 
of  rotation  and  the  difference  in  the  temperature  between  the 
coolant  and  the  wall  (AT)  .  By  increasing  AT  the  heat  transfer 
either  remained  constant  or  increased.  This  increase  was 
greater  for  walls  with  turbulators,  although  the  percentage 
was  not.  This  was  a  result  of  the  higher  heat  transfer  to 
walls  with  turbulators.  When  these  increases,  which  could  be 
larger,  were  compared  to  their  larger  initial  values,  smaller 
percentage  increased  were  found.  Increasing  the  Rossby  number 
induced  higher  heat  transfer  on  the  trailing  surface  while  at 


13 


the  same  time  decreasing  the  leading  surface  value.  They  also 
developed  the  buoyant  parameter  (BP)  such: 

BP-iiE)  (2b)  (BSf)  (3) 

p  V  V 

where  Ap  is  the  pressure  differential,  p  is  the  density  of  the 
coolant,  £2  is  the  angular  velocity,  W  is  the  eccentricity  of 
the  mid  cross-sectional  plane  of  the  duct,  V  is  the  mean 
coolant  velocity  along  the  duct  and  Dh  is  the  hydraulic 
diameter  of  the  duct.  By  increasing  this  value,  a  higher  heat 
transfer  occurred  on  the  trailing  surface,  while  the  leading 
surface  initially  experienced  a  decrease  with  a  subsequent 
increase  when  BP  approached  the  0.2  to  0.4  range.  In  general 
they  found  that  the  locus  of  maximum  heat  transfer  ratios  for 
each  Rossby  number  was  a  function  of  BP. 

Taslim,  Rahman  and  Spring  [Ref  19]  also  confirmed  that  an 
increase  of  heat  transfer  occurred  for  both  the  rotating  and 
stationary  cases  in  a  single  pass  channel  when  staggered  rows 
of  turbulators  were  added  on  opposite  walls.  The  main  thrust 
of  this  experiment  was  to  explore  the  effect  of  Coriolis 
forces  on  internal  heat  transfer.  In  addition,  they 
investigated  the  effect  of  varying  the  blockage  ratio,  e/Dh, 
where  e  denotes  the  turbulator  height  and  Dh  designates  the 
hydraulic  diameter  of  the  passage.  Their  experiments  found  a 
blockage  ratio  of  0.1333  produced  a  maximum  heat  transfer 
coefficient  increase  of  45%  and  a  ratio  of  .333  generated  a 
minimum,  a  6%  decrease.  Guidez  [Ref  20]  used  a  rotating 
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channel  to  correlate  Reynolds  number  and  Nusselt  number  to 
heat  transfer  coefficient.  His  results  indicated  that  the 
pressure  (or  trailing)  side  of  the  channel  experienced 
slightly  increased,  heat  transfer  as  the  Reynolds  number 
increased  while  the  suction  side  heat  transfer  was  weaker.  He 
determined  the  turbulence  at  the  higher  Reynolds  number  became 
larger  than  the  turbulence  induced  by  the  rotation.  His  ratio 
of  rotating  Nusselt  number  to  the  stationary  Nusselt  number 
increased  from  1.0  to  1.75  as  the  channel  accelerated  from  0 
to  600  radians/sec  on  the  pressure  side.  The  ratio  on  the 
suction  side  decreased  from  1.0  to  0.7  and  reached  a  steady 
state  value  at  approximately  175  radians/sec. 

Wagner  et .  al.  [Ref  21]  performed  an  in-depth  experiment 
to  determine  the  effects  that  forced  convection,  Coriolis, 
buoyancy  and  flow  direction  had  on  heat  transfer.  For  outward 
flow,  changes  in  density  ratio  and  Rossby  number  significantly 
affected  the  heat  transfer.  Altering  the  buoyancy  parameter 
for  this  flow  caused  large  changes  in  the  heat  transfer  for 
the  high  pressure  side,  while  changes  made  during  inward  flow 
produced  minimal  changes  in  the  heat  transfer.  Density  ratio 
and  Rossby  number  also  contributed  minimally  during  inward 
flow.  The  suction  side  of  the  channel  proved  relatively 
immune  to  changes  in  heat  transfer  caused  by  modifying  the 
density  ratio,  Rossby  number  or  the  flow  direction.  The  heat 
transfer  was  more  a  function  of  the  buoyancy  than  any  other 
parameter . 
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Morris  and  Hanami-Nasr  [Ref  22]  found  similar  results  in 
the  actions  of  outwardly  directed  flow.  They  found  enhanced 
heat  transfer  on  the  pressure  wall  resulted  from  Coriolis 
induced  secondary  flow  and  a  reduction  on  the  suction  face. 
They  determined  that  turbine  blades  designed  from  heat 
transfer  measurements  made  on  stationary  models  could 
experience  failures  attributable  to  the  inability  to  predict 
local  hot  spots. 

Investigating  the  effects  of  channel  and  turbulator 
geometry  on  internal  heat  transfer,  Taslim  et .  al.  [Ref  23  and 
24]  investigated  a  spanwise  rotating  channel.  In  both  of 
these  experiments,  two  walls  of  a  channel  contained 
turbulators  arranged  at  an  angle  of  attack  to  the  direction  of 
flow  (a)  of  45°  in  a  criss-cross  arrangement.  In  both 
experiments,  blockage  ratios,  aspect  ratios  and  Reynolds 
numbers  were  varied  in  an  attempt  to  determine  their  effect  on 
heat  transfer  coefficients.  Using  the  Rotation  number,  which 
is  defined  as; 


where  Q  is  the  angular  velocity,  Dh  is  the  hydraulic  diameter 
and  Um  denotes  the  mean  velocity  of  air,  the  results  of  both 
investigations  showed  rotational  effects  to  be  more  pronounced 
when  the  test  sections  possessed  higher  aspect  and  lower 
blockage  ratios.  Specifically,  a  steady  increase  in  heat 


transfer  coefficient  observed  on  the  trailing  side  as  the 
Rotation  number  increased  while  the  heat  transfer  coefficient 
on  the  leading  side  showed  a  steady  decrease. 

C .  COMPUTATIONAL  INVESTIGATIONS 

Taylor  et .  al .  [Ref  25],  Moore  and  Moore  [Ref  26],  Medwell 
et .  al .  [Ref  27],  Patankar  and  Prakash  [Ref  28],  Rowley  and 
Patankar  [Ref  29]  and  Patankar  and  Kelkar  [Ref  30]  all 
reported  results  relevant  to  the  next  phase  of  the  present 
investigation:  a  computational  fluid  dynamic  analysis  of 
cooling  airflow  inside  a  turbine  blade.  All  used  a  finite 
volume  analysis  which  produced  good  correlation  to  either 
problems  they  ran  concurrently  or  to  rudimentary  experiments 
performed  previously.  In  addition,  by  including  the  conduction 
in  the  solid  boundary  and  avoiding  upwind  techniques,  they 
concluded  that  this  procedure  will  prove  to  be  exceedingly 
successful  with  more  complex  shapes  in  the  future. 

D.  PURPOSE 

Numerous  experiments  have  been  performed  to  investigate 
the  flowfield  associated  with  an  internally  cooled  turbine 
blade.  There  have  been  considerably  fewer  attempts  to  analyze 
this  flow  using  numerical  techniques.  The  purpose  of  the 
present  analysis  is  to  use  Phoenics,  a  PC-based  computational 
fluid  dynamics  program,  to  assist  in  predicting  and 
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understanding  the  flowf : eld  based  upon  the  experimentations  of 
Taslim  and  Spring  [Ref  7] . 
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II.  THEORY 


A.  GOVERNING  EQUATIONS 

For  a  Newtonian  fluid  with  no  body  forces,  the  law  of 
Conservation  of  Mass  becomes  the  continuity  equation,  Newton's 
Second  Law  yields  the  momentum  equation  and  the  First  Law  of 
Thermodynamics  produces  the  energy  equation.  To  correctly  and 
fully  describe  fluid  flow  it  is  often  necessary  to  include 
additional  equations  designed  to  model  turbulence  in  the 
analysis.  Finally  additional  equations  are  necessary  to 
characterize  the  relationship  between  various  fluid 
properties.  An  example  of  such  an  equation  is  the  equation  of 
state . 

1.  The  Continuity  Equation 

When  considering  fluid  flow  through  a  fixed  control 
volume,  the  continuity  equation  is: 

|fi+V-(pV)=0  (5) 

where  p  is  the  fluid  density  and  V  with  the  overhead  arrow  is 
the  total  velocity  vector  and  u  and  v  are  velocity  components 
in  the  x  and  y  directions  for  two  dimensional  flow.  Using  the 
substantial  derivative 
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(6) 


equation  5  becomes 


■DO  .30 
Dt  dt 


+  V-V() 


^+p(^-v)  =0 


(7) 


For  incompressible  flow 


£E=o 

Dt 


(8) 


and  equation  7  for  incompressible  two-dimensional  flow  becomes 


dx  ay 


(9) 


2.  The  Momentum  Equation 

The  momentum  equation  for  two-dimensional  flow  and  in 
conservative  form  is  [Ref  31] 


dpu  +  dpu2  dpuv-+  dp  _  _  ^xy  _Q 

dt  dx  dy  dx  dx  dy 


(10) 


dpv^  dpuv^  dpv2  __  dp  _&txy  _fayy_n  (11) 

dt  dx  dy  dy  dx  dy 

Here  the  viscous  effects  of  the  flow  are  described  using  the 
shear  stress  terms  Txx,  and  where 
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(12) 


v»»|)  ,19) 

apply.  In  equations  10  through  19,  u  and  v  are  the  velocity 
components  in  the  x  and  y  direction  respectively  ,  p  is  the 
fluid  density,  p  is  the  fluid  pressure  and  (J.  is  the  fluid 
viscosity.  All  terms  with  an  over  bar  indicate  a  time 
averaged  quantity,  while  terms  with  a  primed  superscript 
indicate  fluctuation  quantities. 


21 


3 .  The  Energy  Equation 

Applying  the  First  Law  of  Thermodynamics  to  a  fluid 
passing  through  a  control  volume  yields 

pcij£«i2e+V(JcVD  (20) 

iy  u  XJ  u 

where  the  total  energy  is  comprised  only  of  internal  and 
kinetic  energy.  The  form  of  equation  20  for  heat  transfer  in 
incompressible  two-dimensional  flows  is  written  in  terms  of 
static  temperature  as 


pc,f^pci,^<u  r>*Pc^<v  n- 


6y 


|P  .u-jji ®  ♦ 

3t  dx  By  Bx  By 


(21) 


-k  <*%-*&*>  (*f -pc,^) 

where  the  dissipation  function,  <l>  overbar,  is 


«=— (iH)  +—  ( +ir)  +—(-^) 

~  3.x  x*^3y  3x'  Tyy^  3y; 


(22) 


and  the  fluctuating  dissipation  function  has  been  ignored. 
[Ref  32 j 

4 .  Turbulence  Modeling 

The  unsteady  Navier-Stokes  equations  are  believed  to 
describe  turbulent  fluid  flow.  Currently  the  primary 
technique  employed  in  computational  fluid  dynamic  predictions 
of  fluid  flow  and  heat  transfer  is  to  use  the  time-averaged 
Navier-Stokes  equations.  These  take  the  randomly  generated 
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flow  variables  and  decompose  them  into  their  time  averages  and 
their  instantaneous  fluctuations.  The  time  averages  in  this 
report  are  indicated  with  the  overbar  symbol  while  the 
fluctuations  are  given  the  apostrophe  superscript.  These  new 
terms  are  interpreted  as  apparent  stress  gradients  and  heat- 
flux  quantities  associated  with  turbulent  motion.  These 
quantities  can  then  be  related  to  the  mean  flow  variables 
through  turbulence  models. 

For  the  present  analysis,  the  turbulence  model  used  is 
the  two  equation  model,  commonly  referred  to  as  the  k-e  model, 
discussed  by  Launder  and  Spalding  [Ref  33]  .  The  value  for  the 
effective  turbulent  viscosity,  jj,T,  is 

C„i£  (23) 

where  equals  0.09,  k  is  the  kinetic  energy  of  turbulence 
and  e  is  the  turbulence  dissipation  rate.  They  derived  two 
transport  equations  to  define  their  turbulence  model.  The 
turbulence  kinetic  energy  transport  equation  is  commonly 
modeled  as 


P 


dk 

dt 


+pV-Vk=V( 


VM 


(24) 


where  Prk  is  the  Prandtl  number  for  the  turbulent  kinetic 
energy  and  is  assumed  to  be  1.0.  The  other  transport  equation 
necessary  to  complete  this  turbulence  characterization  is  the 
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one  for  the  turbulence  dissipation  rate.  This  purely 
empirical  formula  is  given  as: 

p|£+py.Ve=V-(-^LVe)+ce  jvf-G-c,  P-^  {25) 

y  dt  y  Prt  **^Tk  k 

In  both  of  these  equations  several  similarities 
appear.  The  left  side  of  both  equations  equate  to  the 
particle  rate  of  increase  of  the  variable.  The  first  term  on 
the  right  side  is  the  diffusion  rate  while  the  second  and 
third  terms  are  the  production  and  dissipation,  respectively. 

This  turbulence  model  proves  more  successful  in  the 
outer  regions  of  the  flow  field  as  opposed  to  the  inner  region 
where  convective  transport  is  relatively  unimportant.  For 
this  region  wall  functions  have  been  developed  for  the  area  in 
the  range  of  30  <  y+  <  200  where  y+  is  defined  as 

1 

2  (26) 

with  xw  being  the  shear  stress,  pw  the  density  and  vw 
representing  the  kinematic  viscosity  all  evaluated  at  the 
wall.  Employing  these  wall  functions  in  the  near  wall  region 
eliminates  the  need  for  a  much  finer  grid  in  this  area. 

B.  THE  GENERAL  DIFFERENTIAL  EQUATION 

Patankar  [Ref  34]  and  Patankar  and  Spalding  [Ref  35] 
recognized  the  relationship  of  the  terms  in  the  model 
transport  equations  and  advocated  a  procedure  of  numerical 
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analysis  of  fluid  flow  which  involves  a  general  differential 
equation.  Their  equation  is 

•Sr  (p4>)  V- < p u<t> )  =V*(rV<|>)  +S  (27) 

ot 

where  <j)  represents  the  dependent  variable  and  T  and  S  refer  to 
the  diffusion  coefficient  and  the  source  term,  respectively. 
These  quantities  are  distinctive  for  the  <(>  being  analyzed. 

The  realization  that  the  applicable  transport  equations  of 
momentum,  the  energy,  turbulence  kinetic  energy  and 
dissipation  and  related  coefficients  can  be  thought  of  as 
particular  cases  of  the  general  <|>  equation  is  an  important 
time-saving  step.  Now  for  the  solution  of  a  fluid  flow  and 
heat  transfer  problem,  only  the  numerical  solution  of  equation 
27  need  be  investigated.  By  constructing  the  applicable 
computer  algorithm  based  on  a  pressure  correction  scheme  [Ref 
34]  in  such  a  manner  as  to  repeatedly  solve  this  same  equation 
with  an  appropriate  definition  for  $  and  the  pertinent  values 
of  r  and  S,  considerable  computational  time  and  memory  savings 
are  realized. 

C.  BOUNDARY  CONDITIONS 

1.  Boundary  Conditions  at  a  Wall 

No-slip  boundary  conditions  are  used  for  the 
velocities  at  a  solid  wall,  where  they  are  all  set  equal  to 
zero.  Wall  temperatures  are  either  explicitly  set  or  the 
temperature  gradient  at  them  is  specified.  When  the  flux  at 
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the  wall  is  set  equal  to  zero  the  wall  assumes  the  properties 
of  an  adiabatic  wall. 

2.  Boundary  Conditions  at  Xn£low 

The  conditions  at  the  inflow  include  all  the  velocity, 
pressure,  temperature  and  turbulence  intensity.  For  this 
investigation  these  values  were  all  specified  at  the  inlet 
condition . 

3.  Boundary  Conditions  at  Outflow 

The  outflow  of  this  problem  required  only  the 
specification  of  the  pressure  to  atmospheric  conditions.  All 
other  conditions  at  the  outflow  were  extrapolated  by  the 
analysis . 
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III.  THE  PHOENICS  COMPUTATIONAL  PROGRAM 


A.  HOW  THE  PHOENICS  COMPUTATIONAL  PROGRAM  WORKS 

Phoenics  is  a  computer  code  which  simulates  fluid-flow, 
heat-transfer , chemical-reaction  and  related  phenomena.  It  is 
capable  of  solving  three-dimensional,  unsteady  and  turbulent 
flow  problems.  Relying  on  the  validity  of  the  assumptions 
used  to  generate  the  input,  it  has  the  ability  to  simulate  the 
effects  of  interaction  between  such  phenomena  and  the 
surrounding  environment .  Phoenics  was  used  throughout  this 
investigation.  For  each  dependent  variable  <|>,  Phoenics 
generates  as  many  difference  equations  as  cells  within  the 
domain.  Thus,  the  total  number  of  equations  Phoenics  must 
solve  is  x  times  y  times  the  number  of  dependent  variables 
present  in  the  domain.  Not  only  are  these  equations  numerous, 
but  they  are  also  very  strongly  coupled.  In  addition  they  are 
also  characteristically  non-linear  due  to  the  functional 
dependence  of  the  coefficients  and  sources  upon  the  auxiliary 
and  dependent  variables.  To  solve  these  non-linear  equations, 
Phoenics  employs  an  iterative  predict-and-correct  technique 
which  reduces  the  imbalance  between  the  left  and  right  sides 
of  each  equation  to  a  negligible  level. 
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B .  GRID  GENERATION 

Phoenics  can  implement  and  use  Cartesian,  cylindrical 
polar  or  curvilinear  or  body-fitted  coordinates.  For  two- 
dimensional  Cartesian  grids,  specification  of  the  length  and 
height  of  the  area  of  analysis  is  required.  The  area  between 
the  x  and  y  boundaries  can  then  be  filled  with  as  many  grid 
lines  as  the  solution  requires.  Spacing  of  these  lines  is 
either  uniform,  exponential  (with  continually  increasing  or 
decreasing  spacing) ,  or  in  a  "power"  arrangement,  which  both 
increases  and  then  decreases  the  spacing  of  the  grid  lines. 
These  last  two  options  allow  the  generation  of  a  grid 
permitting  clustering  of  the  grid  points  closer  together  near 
the  boundaries.  When  different  regions  are  designated,  the 
spacing  can  vary  from  region  to  region  as  necessary.  In 
addition,  by  stipulating  various  regions  in  the  x  or  y 
direction,  irregular  shapes  can  be  analyzed  by  placement  of 
zero-porosity  obstructions  which  eliminate  the  zones  from  the 
flow  field. 
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IV.  BACKGROUND  ON  THE  TURBULATOR  DUCT  PROBLEM 


A.  EXPERIENCE  ON  CASTING  PROBLEMS  IN  TURBINE  BLADES 

As  shown  in  Chapter  I,  internally  routed  cooling  air  flow 
pboves  quite  effective  in  permitting  higher  internal 
operating  temperatures  for  gas  turbine  engines  particularly  if 
turbulators  are  used  in  the  coolant  passages.  The  size,  shape 
and  orientation  of  the  duct  and  the  geometry  of  the 
turbulators  within  these  ducts  is  crucial  to  the  convective 
removal  of  heat  from  the  blades.  As  many  investigations 
showed,  altering  the  aspect  ratio  of  the  duct  or  the 
turbulators,  the  spacing  of  the  turbulators,  their 
orientation,  the  turbulator  pitch-to-height  ratio  and  their 
shape,  greatly  impacts  their  effectiveness  in  promoting  heat 
transfer . 

With  the  continual  emphasis  on  reduction  of  weight  for 
today's  aircraft  and  their  components,  engine  dimensions,  both 
externally  and  internally,  are  also  decreasing.  For  this 
reason  the  casting  process  by  which  turbine  blades  are 
manufactured  is  placing  restrictions  on  the  size  and  shape  of 
the  turbulators  which  can  be  employed.  During  the  turbine 
blade  casting  process  one  of  the  problems  encountered  involves 
the  inability  of  the  molten  materials  to  totally  fill  the 
grooves  forming  the  turbulator. 
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A  second  problem  is  the  tolerances  provided  by  the 


manufacturing  process.  Current  turbine  sizes  require 
turbulators  measuring  approximately  0.010"  in  height.  The 
lower  limit  for  associated  dimensional  tolerances  is  0.002"  - 
0.003".  Thus  a  turbulator  designed  for  a  height  of  0.010"  can 
have  a  variance  of  up  to  30%.  This  difference  can 
significantly  effect  the  heat  transfer  coefficient  for  an 
engine  where  a  0.010"  turbulator  presents  a  high  blockage 
ratio . 

B.  AN  EXPERIMENTAL  INVESTIGATION  INTO  A  TURBULATED  DUCT 
1.  Background  and  Duct  Construction 

Taslim  and  Spring  [Ref  7]  performed  a  series  of 
experiments  to  determine  the  variation  in  heat  transfer  caused 
by  the  geometrical  changes  characteristic  of  small  cast  air¬ 
foils.  Their  method  involved  using  a  46  inch  long  channel 
constructed  of  plexiglas  on  three  sides,  with  the  fourth  side 
made  of  a  polyurethane  slab  incorporating  a  series  of  heaters 
and  liquid  crystal  sheets.  The  width  of  the  heated  wall  was 
3",  which  fixed  that  dimension  for  the  duct.  The  length  of 
the  side  walls  were  changed  to  alter  the  aspect  ratio  during 
the  experiment.  Experiments  were  run  with  aspect  ratios  for 
the  duct  alternating  between  0.50  and  0.55. 

The  first  14  inches  of  the  duct  had  no  turbulators  and 
simulated  the  unturbulated  area  at  the  root  of  the  airfoil, 
while  the  final  32  inches  contained  a  staggered  array  of 
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turbulators  on  the  heated  face  and  the  opposing  plexiglass 
wall.  During  the  experiment,  turbulator  height,  e,  center-to- 
center  pitch,  S,  and  turbulator  width,  w,  were  varied. 
Blockage  ratios,  e/Dh  where  Dh  is  the  hydraulic  diameter,  ran 
between  0.15  and  0.285,  while  pitch-to-height  ratios.  S/e, 
varied  between  4.27  and  14.2. 

2.  Conduct  of  the  Experiment 

After  calibrating  the  liquid  crystal  sheets  in  a  water 
bath  so  a  given  temperature  produced  a  known  color,  and  with 
the  desired  turbulator  and  duct  configuration  established,  air 
flow  was  introduced  into  the  duct  with  a  set  mass  flow  rate. 
This  combination  of  duct  geometry  and  mass  flow  rate 
established  the  Reynolds  number  for  the  test  run. 

Selecting  an  area  approximated  by  x/Dh  =  16,  where  x 
is  the  distance  from  the  test  section  entrance  and  Dh  is  the 
hydraulic  diameter  of  the  unturbulated  section,  the  heat  flux 
was  varied  by  increasing  or  decreasing  the  power  on  the  wall 
heaters  until  a  reference  color  was  obtained  on  the  downstream 
side  of  a  turbulator  near  the  selected  area.  A  picture  of 
this  area  was  taken  and  the  power  altered  until  the  same 
reference  color  was  achieved  in  the  next  adjacent  area.  This 
process  was  repeated  until  the  entire  area  between  the  pair  of 
turbulators  was  covered. 

The  final  experimental  step  required  digitizing  each 
picture  and  measuring  the  area  covered  by  the  reference  color. 
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This  allowed  production  of  an  area-weighted  average  heat 
transfer  coefficient. 

3 .  Results 

The  first  result  desired  from  the  experiment  was  to 
determine  the  impact  of  alleviating  the  problems  caused  by  the 
aforementioned  non-fill  difficulty  by  reducing  the  aspect 
ratio  of  the  turbulator,  ARt .  Taslim  and  Spring  [Ref  7]  found 
a  reduction  of  AR,.  from  two  to  one  reduced  the  Nusselt  number 
on  the  turbulated  wall  by  33%  for  high  blockage  ratio 
turbulators  (e/Dh  =0.22)  with  a  spacing  of  S/e  of  5 .  As 
expected,  a  reattachment  point  formed  directly  downstream  of 
the  turbulators.  This  was  the  point  of  maximum  heat  transfer 
coefficient,  which  then  decreased  monotonically  in  the 
downstream  direction  until  it  approached  the  next  turbulator, 
where  it  starced  to  increase  again.  The  conclusion  from  this 
analysis  was  that  an  optimum  spacing  between  turbulators 
existed  which  allowed  sufficient  length  for  the  flow  to 
reattach,  and  then  a  minimum  of  spacing  thereafter  to  the  next 
turbulator.  The  conclusion  was  that,  for  a  square  turbulator 
(ARt  =  1)  with  a  relatively  high  blockage  ratio,  the  optimum 
S/e  spacing  was  approximately  8.5. 

After  this  determination,  the  effects  of  rounding  the 
turbulators  at  the  tops  and  decreasing  the  sharpness  of  the 
fillets  were  investigated.  Testing  was  performed  on  four  sets 
of  turbulators.  Two  sets  were  manufactured,  both  with  an  ARt 
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of  1.0,  with  one  set  having  sharp  fillets  and  sharp  corners 
and  the  other  having  rounded  fillets  and  rounded  corners.  Two 
similar  sets  were  constructed  with  an  ARC  of  2.0,  with  the 
sharp  or  rounded  fillets  and  corners.  As  expected  the 
rounding  decreased  the  Nusselt  numbers  on  the  turbulated  walls 
for  both  AF-'s,  but  the  ARt  =  2.0  conf iguration  produced  a  17% 
reduction  over  the  sharp  cornered  turbulators  with  the  same 
ARt,  while  the  difference  for  the  turbulators  with  an  ARj  of 
1.0  was  only  5%.  Concurrently  incorporated  was  a  reduction  of 
8.6%  of  surface  area  for  the  turbulator  with  an  ARC  of  2  and 
a  7.15%  reduction  for  the  other  set.  The  contrast  was 
attributed  to  the  difference  in  spacing  between  the  sidewalls 
of  the  adjacent  turbulators.  Although  the  S/e  was  held 
constant,  the  lower  ARC  turbulators  were  wider  and  hence  less 
distance  existed  between  sidewalls.  Thus  for  higher  aspect 
ratio  turbulators,  the  flow  had  sufficient  length  to  reattach 
after  the  turbulator,  while  for  the  lower  aspect  ratio, 
reattachment  did  not  occur.  In  the  latter  case  (tip-to-tip 
turbulator  flow)  there  was  less  dependence  on  turbulator 
geometry . 

From  these  results  an  optimum  turbulator  design  could 
be  attempted.  For  manufacturability,  the  corners  would  be 
rounded  and  the  base  would  be  wider.  For  maximum  heat 
transfer  the  tip  would  be  kept  narrow.  The  resulting  rounded 
trapezoidal  shape  is  termed  a  "jersey-barrier".  At  S/e  of 
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five  the  jersey-barrier  shaped  turbulators  exhibited  the  same 
Nusselt  number  as  high  ARt  turbulators. 

Another  area  investigated  was  the  impact  of 
manufacturing  tolerances  on  the  heat  transfer  characteristics 
of  turbulated  ducts.  Using  a  jersey-barrier  shaped  turbulator 
with  a  nominal  height  of  0.010“  and  a  tolerance  of  ±  0.003"  on 
S/e  of  10,  a  ±40%  deviation  in  the  turbulated  wall  Nusselt 
number  was  experienced  when  dimensions  were  allowed  to 
approach  their  minimum  and  maximum  tolerance  values.  Reducing 
the  S/e  ratio  to  5,  caused  a  variation  from  the  nominal  case 
of  +22%  and  -13%  from  the  maximum  to  minimum  profiles. 

The  final  area  examined  was  the  interaction  between 
the  blockage  ratio,  the  pitch-to-height  ratio  and  the 
turbulator  aspect  ratio.  It  was  found  that  the  decrease  in 
Nusselt  number  associated  with  a  reduction  in  e/Dh  could  be 
completely  offset  by  increasing  the  ARt  and  S/e.  In  addition, 
with  a  constant  S/e  close  to  the  optimum  level  and  a  moderate 
blockage  ratio,  decreasing  the  ARt  from  1.0  to  0.625  resulted 
in  only  a  7%  drop  in  the  Nusselt  number.  This  was 
characteristic  for  geometries  with  the  S/e  closer  to  the 
optimum  value  (8.5)  and  a  lower  value  for  the  blockage  ratio. 
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V.  PHOENICS  SIMULATION  OF  THE  TURBULATOR  DUCT  EXPERIMENT 

The  Phoenics  program  was  used  to  compute  the  results  of 
the  previous  turbulated  duct  experiment.  The  analysis  was 
divided  into  two  basic  steps.  The  first  was  to  generate  a 
rectangular  straight  duct  with  one  wall  heated.  By 
controlling  the  Reynolds  number,  through  varying  inflow 
velocity,  the  Dittus-Boelter  correlation  of  NUS= .  0243Re°-8Pr°-4 
was  verified  [Ref  31].  The  simulation  was  two-dimensional. 

The  second  phase  was  the  generation  of  a  turbulated  duct 
to  reproduce  experimentally  obtained  results  through  numerical 
analysis.  Again  a  rectangular  two-dimensional  straight  duct 
was  used  with  a  staggered  arrangement  of  turbulators  on  the 
top  and  bottom  wall.  The  bottom  wall  was  heated. 

A.  DUCT  AND  FLOWFIELD  SETUP 

1 .  Grid  Generation 

a.  Unturbulated  Square  Duct 

The  first  duct  to  be  modeled  was  30.0m  long.  The 
height  of  the  duct  was  1.0m  and  it  contained  no  turbulators. 
The  bottom  wall  was  heated. 

In  the  x  direction  the  entire  duct  was  considered 
one  region.  It  was  divided  into  100  grid  points  which  were 
stretched  in  the  x  direction  using  a  factor  of  1.10.  This 
configuration  clustered  the  grids  near  the  inlet.  By  doing 
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so,  the  grid  cells  were  square  at  the  entrance  to  the  region, 
which  insured  that  the  flow  achieved  its  fully  developed 
profile  prior  to  entering  the  zone  where  it  was  analyzed. 

In  the  y  direction  there  was  only  one  region. 
There  were  50  grid  points  along  the  y  axis  and  they  were 
clustered  along  the  top  and  bottom  walls. 
b.  Turbulated  Rectangular  Duct 

In  the  second  phase  two  ducts  were  generated.  The 
first  duct  was  1.17m  long  with  the  initial  0.356m 
unturbulated.  The  turbulators  had  an  aspect  ratio,  ARt ,  of 
1.0,  with  a  height,  e,  of  O.C127m  and  a  turbulator  center-to- 
center  pitch,  S,  of  0.0635m.  All  the  turbulators  were  square, 
with  sharp  corners  and  fillets.  The  first  turbulator  was  on 
the  top  wall  and  all  were  positioned  90°  to  the  flow.  The 
second  duct  was  0.654m  long,  with  the  same  unturbulated  region 
at  the  inlet.  The  duct  aspect  ratio  remained  the  same.  Ten 
turbulators  of  the  same  ARt,  S,  and  e  were  again  positioned 
90°  to  the  flow,  with  the  first  one  on  the  top  wall.  The 
bottom  wall  of  both  ducts  in  the  second  phase  was  heated. 

The  longer  duct  was  divided  into  30  regions  along 
the  x  axis.  It  included  175  grid  points  which  were  stretched 
by  a  factor  of  1.2.  Region  two  and  every  fourth  region  there¬ 
after  (6,10,14,...)  through  number  26  were  the  regions  designa¬ 
ted  for  the  turbulators  on  the  top  wall.  Similarly  region 
four  and  every  fourth  region  following  through  number  28  were 
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the  regions  demarcating  the  locations  of  the  bottom 
turbulators.  All  of  the  regions  designating  the  turbulators 
were  0.0127m  long  and  contained  seven  grid  points.  The  odd 
numbered  regions  through  number  29  specified  the  spaces 
between  the  upper  and  lower  turbulators .  Each  was  composed  of 
nine  grid  points  and  was  0.0191m  in  length.  The  final  region 
was  0.369m  long  and  modeled  outflow  from  the  duct.  No 
turbulators  were  located  here  and  it  contained  101  grid 
points . 

In  the  vertical  direction  three  regions  were 
specified.  The  first  region  was  from  the  bottom  wall  up 
0.0127m,  which  corresponded  to  the  height  of  the  turbulators. 
It  contained  15  grid  points.  The  second  region  extended  from 
the  top  of  the  first  up  0.0165m.  Twenty-five  grid  points  were 
allocated  in  this  region.  The  final  region  in  the  y  direction 
simulated  the  top  wall  turbulators.  Like  region  one  it  was 
0.0127m  tall  and  held  15  grid  points.  Figure  1  shows  a 
schematic  of  the  longer  duct  with  an  enlarged  segment  which 
depicts  the  grid  at  the  end  of  the  first  region  and  the  first 
two  turbulators  in  the  turbulated  region. 

This  grid  most  closely  resembled  the  experiment, 
but  did  not  give  adequate  grid  resolution.  An  additional  duct 
was  generated,  which  had  a  finer  mesh  size.  This  duct  was  a 
shortened  version  of  the  fully  turbulated  one.  As  will  be 
explained,  investigation  revealed  that  the  shorter  length  did 
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not  significantly  alter  the  flowfield  prediction  and  did  allow 
a  greater  concentration  of  grid  points  throughout  the  duct. 
This  final  duct  also  included  an  unturbulated  inlet  region, 
whose  dimensions  and  grid  points  were  the  same  as  the  previous 
duct.  The  shortened  turbulated  region  measured  0.298m  long 
and  consisted  of  ten  turbulators  arranged  in  a  staggered  array 
again  with  the  first  one  on  the  top  wall.  Each  turbulator  has 
fifteen  grid  points  from  front  to  back  and  there  are  20  grid 
points  in  the  space  between  each  turbulator.  This  duct  was 
investigated  using  a  uniformly  spaced  grid  and  stretching 
factors  of  1.05,  1.10  and  1.20  all  of  which  clustered  the 

grids  near  the  front  and  back  of  the  turbulators.  There  were 
325  grid  points  in  this  region  and  a  total  of  500  in  the 
entire  horizontal  or  streamwise  direction. 

The  vertical  direction  was  constructed  the  same 
way  as  the  previous  duct  with  three  regions  of  fifteen, 
twenty-five  and  fifteen  grids  each.  Grid  spacing  here  also 
employed  both  stretched  and  uniformly  spaced  points.  Figure 
2  shows  the  same  area  for  this  duct  as  in  Figure  1.  Figure  2 
uses  the  1.2  level  of  grid  spacing  on  the  increased  level  of 
grid  points , 

2 .  Boundary  Conditions 

Analyses  were  run  at  a  number  of  different  Reynolds 
numbers  and  turbulence  levels.  The  Reynolds  number  is  defined 
as 
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Figure  1  -  Schematic  of  Long  Duct  With  Enlargement  of  the 
First  Two  Turbulators  Placed  in  a  Uniformly  Spaced  Grid 
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Figure  2  -  Grid  Using  a  1.20  Stretching  Factor 
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Re= 


(28) 


puDh 

where  Dh  is  the  hydraulic  diameter.  For  a  rectangular  duct, 
Dh  is  defined  as  4A/P  where  A  is  the  area  of  the  duct  and  P  is 
the  perimeter.  For  a  two  dimensional  analysis  a  length  scale 
is  used.  For  the  purposes  of  this  investigation  the  Dh  used 
was  the  value  for  the  actual  Dh  of  the  duct  as  configured  in 
Taslim's  and  Spring's  experiment.  Using  the  values  for 
density  of  1.161  kg/m3  and  for  viscosity  of  1.589  X  10'5  kg/m- 
s,  Table  1  gives  the  inlet  velocities  in  the  horizontal 
direction  associated  with  the  Reynolds  numbers  which  were 
investigated. 

A  value  of  2%  inlet  turbulence  was  specified  for  a 
Reynolds  number  of  10000.  The  turbulent  kinetic  energy  is 
calculated  as  follows 

k=(^.Tu)x[u]2  (29) 

where  Tu  is  the  turbulence  intensity  and  u  is  the  inlet 
velocity,  the  values  for  the  appropriate  turbulence  levels  to 
keep  k  approximately  equal  at  the  inlet  for  all  Reynolds 
numbers  is  also  given  in  Table  1. 
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REYNOLDS  NUMBERS ,  VELOCITIES  AND  TURBULENCE  LEVELS 


INVESTIGATED 


Reynolds 

Velocity 

Turbulence 

Number 

(m/s) 

(%) 

10000 

2.531 

2.000 

15000 

3.796 

0.889 

20000 

5.062 

0.500 

Table  1 

The  inlet  flow  temperature  was  held  at  300°  K  for  all 
runs .  The  bottom  wall  was  heated  with  a  constant  temperature 
of  400°  K,  and  the  top  wall  was  specified  as  adiabatic.  At 
the  outlet  the  external  static  pressure  was  specified. 

The  turbulators  were  obstructions  in  the  flow  and  so 
they  were  specified  as  non-conducting,  adiabatic  obstructions. 
In  addition  they  were  given  the  characteristics  of  possessing 
friction  and  having  zero-porosity.  No  gravity  forces,  heat 
sources  or  internal  plates  were  simulated  in  this  analysis. 
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VI .  COMPUTATIONAL  RESULTS 


A.  UNTURBULATED  PLANAR  DUCT 

As  an  initial  test  case  an  unturbulated  square  duct  was 
evaluated  by  the  Phoenics  computational  program.  For  fully 
developed  laminar  flow,  a  parabolic  velocity  profile  develops 
with  a  maximum  velocity,  umax=2umean .  For  turbulent  flow  the 
velocity  profile  is  much  flatter  in  the  center  due  to  the 
mixing  motion  and  momentum  exchange  and  u^^l . 2uraean . 

The  values  of  The  Reynolds  numbers  for  this  investigation 
were  10000  and  20000.  The  geometry  of  the  chosen  duct 
provided  a  30:1  length  to  height  ratio.  This  allowed  a 
sufficient  length  for  the  flow  to  develop  fully  prior  to  being 
investigated.  The  zone  chosen  for  the  analysis  of  the  flow 
was  located  between  20  and  25  diameters  downstream  of  the 
inlet . 

Figure  3  is  a  profile  of  the  flow  located  25  diameters 
from  the  inlet.  As  can  be  seen  it  presents  the  classic  shape 
associated  with  turbulent  fluid  flow  in  a  duct.  The  maximum 
value  of  the  velocity  of  the  flow  occurs  at  mid-duct  and  is  a 
value  of  0.1181  m/s.  Using  the  Phoenics  results  file  the  mean 
average  for  velocity  at  this  point  in  the  flow  was  0.998  m/s. 
The  ratio  is  approximately  1.18,  which  is  indicative  of 
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Figure  3  -  Flow  Profile  25  Diameters  Downstream  in  a  Square 
Unturbulated  Duct 


turbulent  flow.  A  plot  of  the  skin  friction  coefficient,  cf, 
where 


(30) 


xv  is  the  shear  stress  along  the  wall  and  S  is  the  area,  on 
the  bottom  wall  (Figure  4)  shows  the  increase  associated  with 
the  viscous  effects  as  the  flow  progresses  along  the  duct. 
The  value  for  this  shear  stress  along  the  top  and  bottom  wall 
for  the  area  of  fully  developed  flow  is  8.175  X  10'5  kPa. 
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Figure  4  -  Skin  Friction  Coefficient  for  the  bottom  Wall  of 
a  Square  Duct 


B.  TURBULATED  RECTANGULAR  DUCT 
1.  Full  Turbulated  Duct 

A  solution  for  the  turbulated  full  duct  was  generated 
first.  As  stated  previously  this  analysis  was  run  with  a 
relatively  coarse  grid  in  order  to  simulate  Taslim's  [Ref  7] 
experiment.  This  run  included  a  Reynolds  number  of  10000  and 
converged  at  approximately  1900  iterations. 

As  shown  in  Figure  5,  at  both  the  first  and  second 
turbulators  the  flow  generated  a  stagnation  point  on  the 
upstream  face  of  the  turbulator  and  recirculation  occurred 
ahead  of  the  turbulator.  Figure  6  is  a  plot  of  the  velocity 
vectors  for  the  area  of  the  duct  adjacent  to  the  second 
turbulator.  The  flow  did  not  recirculate  upstream  of 
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Figure  5  -  x-Direction  Velocity  Vectors  at  Turbulators  1  and 
2 
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Figure  6  -  Velocity  Vectors  Between  Turbulators  2  and  4 
Using  a  Constant  Grid  Spacing 
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turbulators  three  and  four,  because  of  the  influence  of  the 
upstream  turbulators  two  and  three  respectively.  This  flow 
pattern  became  repeatable  from  this  location  onward  until  the 
end  of  the  turbulated  section. 

The  position  of  the  reattachment  point  downstream  of 
the  turbulator  was  dependent  on  the  location  within  the  duct. 
As  the  flow  progressed  down  the  duct,  the  reattachment  point 
moved  closer  to  the  downstream  face  of  the  turbulator 
eventually  reaching  a  fixed  distance  downstream  of  the 
turbulator.  The  area  between  turbulators  eight  and  ten  was 
the  first  to  approximate  this  "periodic"  flowfield.  Figure  7 
shows  the  vectors  as  they  appeared  for  the  flow  in  this  area 
and  Figure  8  plots  the  x-direction  velocity  components  for  the 
first  grid  point  away  from  the  lower  surface. 

Another  area  of  interest  was  the  behavior  of  the 
freestream  kinetic  energy  (k)  of  the  turbulence.  The  mid-duct 
value  for  turbulent  kinetic  energy  was  plotted  for  the  entire 
length.  Figure  9  shows  an  initial  free  stream  decay  of 
turbulence  in  the  unturbulated  section,  after  which  the 
turbulence  increased  dramatically  at  the  first  turbulator  and 
experienced  additional  increases  at  each  subsequent  one. 
Initially  each  gain  was  significant,  but  eventually  they 
became  smaller  in  magnitude.  After  the  eighth  turbulator  a 
near  constant  value  was  approximated  for  the  free-stream 
turbulence.  This  asymptotic  value  was  two  orders  of  magnitude 
over  the  minimum  values  experienced  in  the  unturbulated 
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Figure  7  -  Velocity  Vectors  Between  Turbulators  8  and  10 
Using  a  Constant  Grid  Spacing 
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Figure  8  -  x-Direction  Velocity  Component  Between 

Turbulators  8  and  10  Using  a  Constant  Grid  Spacing 


section.  After  the  last  turbulator,  the  flow  entered  the 
unturbulated  outlet  portion  of  the  duct  and  turbulence  once 
again  decreased  due  to  viscous  dissipation. 


Figure  9  -  Turbulent  Kinetic  Energy  for  a  Full  Turbulated 
Duct 
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Using  the  above  data,  the  tenth  turbulator  was 
employed  as  the  cutoff  point  for  the  full  investigation. 

2.  Abbreviated  Turbulated  Duct 

a.  Uniformly  Spaced  Grid  Simulation 

Using  a  uniformly  spaced  grid  and  a  duct  shortened 
to  ten  turbulators,  a  simulation  was  run  at  a  Reynolds  number 
of  10000  with  a  turbulence  intensity  of  2.0%.  The  development 
of  the  flow  was  consistent  with  the  full  length  duct.  Both 
the  first  turbulator  (on  the  top  wall)  and  the  second  one  (on 
the  bottom  wall)  developed  regions  of  recirculation  prior  to 
them.  Figure  10  is  a  plot  of  the  velocity  vectors  for  the 
area  prior  to  turbulator  two  and  it  shows  the  flow  pattern 
with  a  distinct  stagnation  point. 

The  flow  for  the  area  downstream  of  turbulator 
two,  however  was  markedly  different  from  the  structure  found 
using  the  full  length  duct.  Here  the  entire  area  between  the 
second  and  fourth  turbulators  exhibited  recirculation  of  the 
flow.  No  reattachment  point  along  the  bottom  wall  existed. 
Figure  11  is  a  plot  of  the  vectors  in  this  area  and  shows  the 
total  reverse  flow.  A  plot  of  the  horizontal  velocity  at  the 
first  grid  line  off  the  bottom  wall  for  this  region  shows  the 
negative  value  of  this  component  in  Figure  12. 

After  the  fourth  turbulator,  the  flow  along  the 
bottom  wall  begins  to  exhibit  the  characteristics  of  the  pre¬ 
vious  duct.  The  flow  over  the  turbulator  turns  down  to  the 
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Figure  10  -  x-Direction  Velocity  Vectors  Prior  to 

Turbulators  1  and  2  Using  a  Constant  Grid  Spacing 
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Figure  11  -  Velocity  Vectors  Between  Turbulators 
Using  a  Constant  Grid  Spacing 
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reattachment  point  establishing  a  recirculation  zone 
downstream  of  the  turbulator.  The  flow  then  turns  toward  the 
next  downstream  turbulator.  On  the  upstream  face  of  this 
turbulator  no  recirculation  occurs  and  the  flow  proceeds  over 
the  top  with  no  separation  on  the  upper  surface. 


Figure  12  -  Horizontal  Velocity  Component  Between 

Turbulators  2  and  4  Using  a  Constant  Grid  Spacing 


Using  data  generated  in  the  Phoenics  results 
files,  the  area  between  the  eighth  and  tenth  turbulators  was 
investigated  further.  Figures  13  and  14  show  the  horizontal 
and  vertical  velocity  components  respectively  of  the  flow  in 
this  area  along  the  first  grid  line  off  the  bottom  wall.  The 
horizontal  velocity  shows  a  minimum  value  of  -0.51  m/s  in  the 
recirculation  zone  with  the  reattachment  point  located  0.014m 
downstream  of  the  turbulator.  The  flow  accelerates  for  0.027m 
where  it  reaches  a  maximum  of  1.9  m/s.  After  the  reattachment 
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point  all  values  for  the  flow  are  positive,  verifying  the  lack 
of  a  recirculation  region  on  the  face  of  the  tenth  turbulator. 

Figure  14  illustrates  the  vertical  velocity 


component's  maximum  values  are  located  in  those  areas  immedi¬ 
ately  after  and  immediately  prior  to  the  turbulator  faces. 


Figure  13  -  Horizontal  Velocity  Component  Between 

Turbulators  8  and  10  Using  a  Constant  Grid  Spacing 


Figure  15  shows  a  plot  of  the  calculated  Stanton  number  versus 
locacion  in  the  duct.  The  first  peak  is  located  adjacent  to 
the  downstream  face  of  turbulator  eight  and  coincides  with  the 
maximum  vertical  component  of  flow  and  a  near  zero  value  of 
horizontal  velocity.  Due  to  the  location  of  this  peak 
adjacent  to  the  wall  of  the  turbulator,  predictions  of  global 
heat  transfer  were  not  consistent  with  those  found  in  the 
experiment.  This  will  be  examined  closer  in  the  next  section 
using  grid  stretching  to  cluster  more  grids  near  the  wall. 
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Figure  14  -  Vertical  Velocity  Component  Between  Turbulators 
8  and  10  Using  a  Constant  Grid  Spacing 


Figure  15  -  Stanton  Number  Between  Turbulators  8  and  10 
Using  a  Constant  Grid  Spacing 

The  next  three  peak  Stanton  numbers  occur  between 
0.010m  and  O.Olirn  from  the  downstream  face  of  the  turbulator. 
When  compared  to  the  velocity  component  plots,  two  distinct 
peaks  are  found  to  straddle  the  reattachment  point  with  one 
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rather  elongated  peak  occurring  over  the  reattachment  zone. 
The  final  peak  occurs  near  the  upstream  face  of  the  tenth 
turbulator.  This  maximum  is  located  at  the  position  of  the 
maximum  vertical  component  of  the  flow  and  where  the 
horizontal  component  is  nearing  zero. 

The  effect  of  the  fluid  flow  and  heat  transfer 
within  the  duct  is  shown  in  Figure  16  which. is  a  plot  of  the 
wall  temperature  between  the  two  turbulators.  The  maximum 
temperature  develops  near  the  downstream  face  of  the  eighth 
turbulator.  The  temperature  first  decreases  and  then  in¬ 
creases  to  a  relative  maximum  prior  to  the  reattachment  point. 


Figure  16  -  Temperature  Between  Turbulators  8  and  10  Using 
a  Constant  Grid  Spacing 


At  the  reattachment  point  a  relative  minimum  temperature  re¬ 
sults  from  the  predominantly  downward  flow.  Another  relative 


maximum  occurs  on  the  downstream  side  of  the  reattachment 
point  followed  by  the  minimum  temperature  region.  This  region 
corresponds  to  the  area  where  the  horizontal  velocity 
component  is,  in  general,  one  order  of  magnitude  greater  than 
the  ver  tical  component .  A  relative  temperature  maximum  arises 
near  the  upstream  face  of  the  tenth  turbulator,  corresponding 
to  the  relative  Stanton  number  maximum  in  the  same  location. 

The  development  of  the  kinetic  energy  of 
turbulence  is  plotted  in  Figure  17.  An  initial  decrease  in 
the  unturbulated  area  of  the  duct  approaches  a  minimum  of  4.3 
X  10'2.  As  the  flow  enters  the  turbulated  region  the  effects 
of  each  turbulator  is  observed.  The  characteristic  increases 
experienced  by  KE  at  each  turbulator  initially  are  on  the 
order  of  3.5  to  4.0.  Eventually  they  decline  to  approximately 
0.1.  The  total  increase  of  KE  approximates  two  orders  of 
magnitude . 

b.  Simulation  Using  Grid  Stretching 

Using  the  same  duct  as  above  but  with  a  1.2  factor 
to  stretch  the  grid  in  all  regions  in  both  the  x  and  y 
direction,  another  simulation  was  run.  Again  the  flow  was  run 
at  a  Reynolds  number  of  10000  with  an  inlet  turbulence 
intensity  of  2.0%.  The  initial  development  of  the  flow 
remained  the  same.  The  characteristic  recirculation  regions 
formed  prior  to  both  turbulators  one  and  two.  The  area 
between  turbulators  two  and  four  again  contained  totally 
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Figure  17  -  Turbulent  Kinetic  Energy  for  the  Shortened  Duct 

separated  flow  with  no  reattachment  point.  Subsequent  flow 
developed  in  a  pattern  similar  to  the  flow  associated  with  a 
uniformly  spaced  grid.  A  reattachment  point  developed  after 
the  turbulator  with  an  area  of  recirculation  upstream  of  it 
between  the  turbulator  and  the  reattachment  point .  No 
recirculation  or.  the  upstream  face  occurred  for  any  flow  after 
the  second  turbulator. 

Again  using  the  Phoenics  results  file  the  area 
between  turbulators  eight  and  ten  can  be  examined  closely. 
Figures  18  and  19  show  the  horizontal  and  vertical  components 
of  the  flow  velocity  respectively  in  the  analyzed  regions. 
Figure  18  indicates  this  grid  generated  a  larger  recirculation 
zone  than  the  uniformly  spaced  grid.  The  stretched  grid 
developed  a  recirculation  zone  of  0.017m  versus  the  uniform 
grid's  0.014m.  The  stretched  grid  also  generated  smaller 


Figure  18  -  Horizontal  Velocity  Component  Between 


Turbulators  8  and  10  Using  a  1.20  Stretching  Factor 


8  and  10  Using  a  1.20  Stretching  Factor 


maximum  and  minimum  values  for  the  vertical  component  of  the 
velocity  although  the  locations  were  relatively  similar. 
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The  plot  of  the  Stanton  number  between  the  two 


turbulators  (Figure  20)  showed  considerable  differences  when 
coupled  with  the  one  for  the  uniform  grid  (Figure  15) . 


Figure  20  -  Stanton  Number  Between  Turbulators  8  and  10 


Using  a  1.20  Stretching  Factor 

The  stretched  grid  solution  located  a  relative 
maximum  adjacent  to  the  downstream  face  where  the  vertical 
component  of  the  velocity  was  a  maximum,  and  the  horizontal 
component  was  near  zero.  It  also  located  two  relative 
maximums  on  either  side  of  the  reattachment  point  as 
determined  from  Figure  16.  The  next  point  identified  is  the 
maximum  for  the  Stanton  number  in  the  area  between  these  two 
turbulators.  It  is  positioned  just  downstream  from  the 
minimum  vertical  velocity  component,  approximately  0.021m  from 
the  downstream  face  of  turbulator  8.  Its  effect  is 
discernable  from  Figure  21,  the  plot  of  the  temperature 
distribution  between  the  two  turbulators. 
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Figure  21  -  Temperature  Between  Turbulators  8  and  9  Using  a 
1.20  Stretching  Factor 


The  maximum  temperature  is  located  precisely  where 
the  maximum  heat  transfer  occurs,  as  indicated  by  the  Stanton 
number.  The  final  relative  peak  for  the  Stanton  number  is 
located  proximate  to  the  upstream  face  of  the  tenth  turbu- 
lator.  As  in  the  previous  solution  it  is  collocated  with  the 
maximum  value  for  the  flow's  vertical  component  of  velocity 
and  where  the  horizontal  component  is  approaching  zero. 

Since  the  maximum  Stanton  number  and  temperature 
occurs  near  the  position  on  the  grid  where  the  clustering  of 
the  grid  points  for  the  upper  turbulator  is  greatest,  it  ap¬ 
pears  that  this  maximum  is  grid  related,  rather  than  a  valid 
depiction  of  the  actual  situation.  Another  variation  of  this 
configuration  was  run  using  a  grid  stretching  factor  of  1.10. 
As  seen  in  Figure  22,  the  peak  appearing  in  Figure  20  is  not 
present  for  this  new  configuration.  The  new  peak  is  located 
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Figure  22  -  Stanton  Number  Between  Turbulators  8  and  10 

Using  a  1.10  Stretching  Factor 


on  the  upstream  side  of  the  attachment  point  and  is  half  the 
magnitude  of  that  experienced  for  a  stretching  of  1.20. 

A  final  attempt  to  examine  this  grid  generated 
anomaly  was  made  using  a  grid  stretching  factor  of  1.05. 
While  this  tended  to  cluster  the  grid  closer  to  the  downstream 
face  of  the  turbulator,  which  appeared  to  rectify  the  spike 
from  the  uniformly  spaced  grid,  it  did  not  create  an  artifi¬ 
cial  spike  in  the  area  adjacent  to  or  in  the  reattachment 
zone.  Figure  23  is  the  Stanton  number  plot  for  the  flow  be¬ 
tween  turbulators  eight  and  ten  for  this  grid  configuration. 
The  most  obvious  difference  between  this  and  the  previous 
cases  is  the  reduction  in  the  magnitudes  of  the  peaks.  All 
the  peaks  in  this  plot  which  are  located  at  the  same  approxi¬ 
mate  location  as  the  peaks  associated  with  the  previous  grids 
are  lower  by  an  order  of  magnitude.  The  peak  adjacent  to  the 
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downstream  face  of  the  turbulator  is  also  approximately  an 
order  of  magnitude  smaller  than  the  peak  associated  with  the 
uniformly  spaced  grid.  This  analysis  confirms  that  the 
solution  for  this  area  is  highly  grid  dependent. 


Figure  23  -  Stanton  Number  Between  Turbulators  8  and  10 


Using  a  1.05  Stretching  Factor 

The  plot  of  the  temperature  for  the  area  between 
turbulators  eight  and  ten  using  a  1.05  stretching  factor 
(Figure  24)  shows  that  the  maximum  temperature  is  associated 
with  the  relative  maximum  Stanton  number  adjacent  to  the 
downstream  face  of  the  turbulator.  The  other  three  relative 
maximums  for  temperature  are  also  associated  with  relative 
maximums  of  the  Stanton  number.  Contrary  to  all  previous 
investigations  a  relative  minimum  of  the  temperature  occurs  at 
a  relative  maximum  Stanton  number. 

Defining  shear  stress  according  to  equation  13  and 
plotting  it  for  this  portion  of  the  duct  results  in  Figure  25. 
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Here  the  plot  shows  that  the  value  of  the  shear  stress  at  the 


point  of  the  relative  maximum  Stanton  number  is  zero. 


Figure  25  -  Shear  Stress  Between  Turbulators  8  and  10  Using 


a  1.05  Stretching  Factor 


This  duct's  turbulent  kinetic  energy  plot  (Figure 
26)  is  very  similar  to  that  for  the  previous  duct.  Original 
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values  are  the  same  with  the  same  initial  decrease 
experienced.  Although  they  both  reached  the  same  level  of 
turbulence,  the  stretched  grid  appears  to  take  longer  to 
develop  this  value.  The  final  increases,  however,  are  still 
more  than  an  order  of  magnitude  smaller  than  the  initial  ones 
and  can  be  considered  to  be  an  asymptotic  condition. 


Figure  26  -  Turbulent  Kinetic  Energy  for  the  Abbreviated 


Duct  Using  a  1.05  Stretching  Factor 

Another  investigation  which  was  done  was  the 
determination  of  the  effects  of  varying  the  turbulence  inten¬ 
sity  at  the  inlet.  Figures  27  and  28  show  plots  of  the  KE  at 
a  mid-duct  location  for  a  Reynolds  number  of  20000  with  the 
initial  turbulence  intensity  being  2.0%  and  0.5%  respectively. 
As  expected  the  starting  and  development  areas  of  the  duct  for 
both  conditions  are  different.  However  once  the  flow  enters 
into  the  turbulated  area,  the  flow  patterns  and  the  levels  of 
turbulence  start  to  converge.  From  this  it  can  be  concluded 
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that  the  level  of  turbulent  kinetic  energy  at  an  asymptotic 
location  in  a  duct  is  unaffected  by  the  level  of  turbulence 


intensity  far  upstream. 


of  2000  and  a  Turbulence  Intensity  of  2% 


Figure  28  -  Turbulent  Kinetic  Energy  for  a  Reynolds  Number 
of  20000  and  a  Turbulence  Intensity  of  0.5% 
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VZ.  COMPARISON  OF  RESULTS  WITH  EXPERIMENT 


A.  UNTURBULATED  SQUARE  DUCT 

The  results  of  the  unturbulated  square  duct  can  be 
compared  with  the  results  of  the  Dittus-Boelter  equation  which 
states 

Nu=0 . 0243Pr°-4i?e°-8  (31) 

where  Nu  is  the  Nusselt  number.  This  equation  is  valid  when 
the  Prandtl  number  is  between  0.7  and  160.0  and  the  Reynolds 
number  is  greater  than  104.  From  Phoenics,  the  steady  state 
Stanton  number  was  calculated  for  Reynolds  numbers  of  10000 
and  20000.  This  is  a  measure  of  the  heat  transfer  occurring 
at  the  bottom  wall  of  the  duct .  This  number  was  then 
converted  to  the  Nusselt  number.  Figure  29  is  a  plot  of  the 
Dittus-Boelter  equation  (dashed  line)  and  the  experimental 
data  between  Reynolds  numbers  of  10000  and  20000. 

The  disagreement  between  the  Phoenics  generated  values  for 
Nusselt  number  and  the  Dittus-Boelter  correlation  is 
considerable.  The  slope  is  opposite  between  the  two  and  the 
magnitude  is  one  full  order  higher  for  the  Phoenics  generated 
data.  Operating  at  the  lower  end  of  the  effective  region  for 
the  Dittus-Boelter  correlation  and  grid  resolution  problems 
are  possible  causes  of  this  discrepancy. 
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Figure  29  -  Nusselt  Number  vs  Reynolds  Number  for  a  Square 
Duct 

B.  TURBULATED  RECTANGULAR  DUCT 
1.  Constant  Grid  Spacing 

The  plot  of  the  Stanton  number  for  the  uniformly 
spaced  grid  at  a  Reynolds  number  of  10000  (Figure  15)  is 
dominated  by  the  peak  near  the  downstream  wall.  This  is 
specific  to  this  test  case  and  does  not  appear  on  the  plots 
for  Reynolds  numbers  of  15000  or  20000.  This  adversely 
affects  the  entire  slope  of  the  numerically  generated  Nusselt 
number  solutions. 

Since  each  grid  is  a  constant  width,  the  average 
Stanton  number  over  the  entire  zone  can  be  resolved  by  deter¬ 
mining  the  value  for  the  Stanton  number  at  each  grid  point  and 
averaging  it.  Plotted  on  log-log  paper  with  the  results  of 
Taslim  and  Spring  [Ref  7]  ,  the  results  of  the  uniformly  spaced 
grid  show  a  negative  slope,  while  that  of  the  experimentally 


generated  data  is  positive.  In  addition,  the  experimental 
data  corresponds  closely  to  the  slope  of  the  Dittus-Boelter 
correlation  line.  A  plot  of  these  two  lines  is  shown  in 
Figure  30  with  the  experimental  data  being  depicted  with  a 
dashed  line. 


Figure  30  -  Nusselt  Number  vs  Reynolds  Number  for  a 


Uniformly  Spaced  Grid 

2.  Grid  Spacing  Using  Stretching 

Three  different  factors  of  stretching  were  used 
throughout  this  analysis:  1.20,  1.10  and  1.05.  Since  the 

grids  are  variably  spaced,  the  average  Stanton  number  requires 
each  individual  Stanton  number  be  multiplied  by  the  width  of 
its  grid.  All  of  these  values  must  then  be  added  together  and 
then  divided  by  the  total  width  to  determine  the  average. 
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The  plots  of  the  Stanton  numbers  of  the  1.20  and  1.10 
stretched  grids  are  dominated  by  exceeding  high  artificial 
peaks  (Figures  20  and  22)  .  When  converted  to  Nusselt  numbers, 
averaged  with  the  other  values  and  plotted  against  Reynolds 
numbers,  these  artificially  high  values  induce  negative  slopes 
on  the  resulting  curves.  This  runs  contrary  to  not  only  the 
results  of  Taslim,  but  also  the  Dittus-Boelter  correlation  and 
the  Colburn  correlation  which  states  Nu=f /2Pr1/3Re  where  f/2  is 
a  dimensionless  shear  stress  term  [Ref  36]  .  All  of  these 
clearly  indicate  an  increase  of  Reynolds  number  should  cause 
an  increase  of  Nusselt  number.  A  plot  of  the  resultant  graph 
for  the  stretching  of  the  grid  by  factors  of  1.10  and  1.20  are 
shown  in  Figures  31  and  32  respectively. 


Figure  31  -  Nusselt  Number  vs  Reynolds  Number  Using  a  Grid 
Stretching  Factor  of  1.10 
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Figure  32  -  Nusselt  Number  vs  Reynolds  Number  Using  a  Grid 
Stretching  Factor  of  1.20 

The  stretching  factor  of  1.05  produced  no  such 
anomalies.  The  values  for  Stanton  number  were  more 
consistent,  with  no  one  peak  exceeding  the  others  by  an  order 
of  magnitude  or  more.  The  plot  of  this  resultant  curve  is  not 
similar  to  the  slope  of  that  generated  by  Taslim  and  Spring. 
As  seen  in  Figure  33,  with  its  minimal  slope,  it  more  closely 
approximates  a  constant  Stanton  number  instead  of  an 
increasing  one  with  increasing  Reynolds  number. 

The  calculation  of  the  Stanton  number  by  the  Phoenics 
algorithm  is  a  blind  operation  which  does  not  reveal  the 
equations  or  values  used.  To  determine  if  the  variance  in  the 
Phoenics  generated  Nusselt  numbers  are  the  result  of  a  grid 
resolution  problem,  the  Nusselt  number  was  hand  calculated 


Figure  33  -  Nusselt  Number  vs  Reynolds  Number  Using  a  Grid 


Stretching  Value  of  1.05 


from  the  values  of  component  variables  output  by  Phoenics  and 
compared  with  the  values  computed  in  the  code. 

First  the  heat  flux,  q0'',  is  defined  as 


where  K  is  the  thermal  conductivity  of  air  at  300°K,  CH  is  a 
constant  in  the  k-e  turbulence  model  equal  to  0.09,  Cp  is  the 
specific  heat  of  air  at  constant  pressure  and  ot  is  the 
turbulent  Prandtl  or  Schmidt  number  [Ref  37] . 

Once  the  value  of  the  heat  flux  is  determined,  a  local 
heat-transfer  coefficient,  hx,  is  specified  as 


73 


(33) 


where  t0  is  the  specified  temperature  at  the  bottom  wall.  The 
t„,  temperature  is  located  at  mid-duct.  Finally  the  local 
Nusselt  number  can  be  generated  from 

Nu=h±-X  (34) 

K 

Here  the  x  value  is  taken  as  the  width  of  the  cell  for  which 
the  temperature,  kinetic  energy  and  dissipation  of  turbulence 
are  being  considered. [Ref  38] 

A  plot  of  the  calculated  Nusselt  number  versus  the 
distance  from  the  inlet  is  shown  in  Figure  34.  The  dashed 
line  plot  is  the  Phoenics  generated  values.  A  very  striking 
resemblance  exists  between  the  two  graphs.  All  but  one  of  the 
maximums  from  the  Phoenics  results  are  identifiable  for  the 
calculated  values .  The  one  value  which  does  not  appear  on  the 
chart  of  the  calculated  values  is  the  spike  at  the  reattach¬ 
ment  point.  The  indication  is  this  spike  is  grid  generated. 

As  a  final  comparison,  the  values  for  the  average 
calculated  Nusselt  numbers  for  Reynolds  numbers  of  10000  and 
20000  were  plotted  versus  these  Reynolds  numbers  in  Figure  35. 
Also  plotted  were  the  values  obtained  from  Phoenics  for  a  1.05 
stretching  factor  solution  (solid  line)  and  the  experimental 
values  (dashed  and  dotted  line)  .  The  slope  for  the  calculated 
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DISTANCE  FROM  INLET  (m) 


Figure  34  -  Calculated  and  Phoenics  Generated  Values  for 
Nusselt  Number  Between  Turbulators  8  and  10 

values  is  much  closer  to  the  experiment  than  that  produced  by 
Phoenics,  but  the  magnitude  is  much  greater. 

Because  of  the  close  correlation  between  the  Phoenics 
generated  program  and  the  hand  calculated  one,  the  failure  to 
produce  Taslim's  results  over  a  range  of  Reynolds  numbers 
occurs  due  to  an  inability  to  adequately  resolve  the  flowfield 
with  this  grid  size. 
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Figure  35  -  Nusselt  Number  vs  Reynolds  Number  Using 

Calculated,  Phoenics  Generated  or  Experimental  Values 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  analysis  of  the  data  presented  in  this  investigation 
leads  to  a  number  of  conclusions  about  airflow  in  a 
rectangular  duct  with  square  turbulators  placed  in  a  staggered 
array.  Using  a  S/e  ratio  of  5,  an  e/Db  of  0.235,  an  ARt  of  1.0 
(Taslim's  value  of  Dh)  ,  with  an  e  of  0.0127m,  the  following 
can  be  concluded: 


•  Flow  approaching  the  first  set  of  turbulators  will 
experience  recirculation  on  the  upstream  side  of  only  the 
first  pair  of  turbulators.  All  turbulators  thereafter 
will  direct  the  flow  up  and  over  the  turbulator  and  will 
not  generate  a  separation  point  for  the  recirculation  of 
the  oncoming  flow. 

•  Flow  over  the  top  of  the  turbulators  will  remain  attached 
and  create  no  recirculation  zones.  This  was  true  for  a 
range  of  Reynolds  numbers  between  7000  and  30000  with 
turbulence  intensity  varying  between  0.22  and  4.1. 

•  No  reattachment  point  will  develop  for  flow  over  the  first 
or  second  turbulator.  The  entire  region  between 
turbulators  two  and  four  will  be  recirculating  flow.  This 
is  the  only  area  of  the  duct  which  will  experience  this 
flow  structure. 

•  The  pattern  of  flow  between  all  sets  of  turbulators  after 
the  initial  two  on  either  wall  will  maintain  the  same 
characteristic  structure.  The  flow  will  develop  a 
recirculation  region  directly  after  the  turbulator 
followed  by  a  reattachment  point.  An  area  of  acceleration 
of  flow  velocity  toward  the  downstream  turbulator  will 
ensue  with  the  flow  traveling  over  the  turbulator.  Each 
succeeding  section  results  in  a  smaller  recirculation  zone 
until  an  asymptotic  dimension  is  reached. 
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•  Maximum  heat  transfer  occurs  at  the  reattachment  point  of 
the  downstream  recirculation  zone.  In  general  four 
relative  maximum  heat  transfer  locations  occur:  on  the 
downstream  face  of  the  turbulator,  on  either  side  of  the 
reattachment  point  and  on  the  upstream  face  of  the 
following  turbulator. 

•  The  kinetic  energy  of  turbulence  throughout  a  turbulated 

duct  maintains  a  distinct  trace.  Each  turbulator 

increases  the  level  as  the  medium  passes  over  it.  Initial 
turbulators  provide  a  significantly  greater  level  of 
increase  than  those  located  further  along  the  duct.  After 
eight  turbulators,  an  asymptotic  level  of  kinetic  energy 
is  reached. 

•  The  asymptotic  level  of  turbulence  in  a  duct  is 

independent  of  the  turbulence  intensity  at  the  inlet.  In 
this  duct,  the  individual  variation  of  kinetic  energy 

caused  by  each  turbulator  increased  or  decreased  with 

respect  to  a  similar  change  in  turbulence  intensity.  Thus 
altering  the  turbulence  intensity  did  not  change  the 
results  of  the  simulation. 

•  The  detailed  results  are  highly  grid  dependent.  The 
choice  of  uniform  or  clustered  grids  not  only  effects  the 
results,  but  the  power  of  the  clustering  can  also  impact 
them.  A  converged  solution  can  exhibit  grid-generated 
phenomena  which  have  no  basis  in  actuality  and  for  which 
the  investigator  must  constantly  inspect. 

•  This  investigation  was  a  2 -dimensional  analysis  of  an 

experiment  performed  in  three  dimensions .  The 

experimental  results  showed  a  variation  in  the  spanwise 
direction  of  the  duct.  The  effects  of  the  walls  on  the 
flow  caused  lower  values  for  the  heat  transfer  coefficient 
near  the  sidewalls  with  an  increase  closer  to  the 
centerline.  A  2 -dimensional  analysis  would  not  have  the 
effects  of  the  sidewalls  included  and  would  produce  this 
centerline  higher  value.  This  possibly  explains  the 
elevated  values  seen  in  the  computed  results. 

B .  RECOMMENDATIONS 

Phoenics  is  a  powerful  program  with  the  ability  to  solve 
highly  complex  flow  problems.  Its  strength  lies  in  the  ease 
with  which  parameter  can  be  altered  and  reevaluated  for  a 
given  problem.  The  present  analysis  was  done  using  strictly 
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a  Cartesian  coordinate  system  in  two  dimensions.  With 
Phoenics'  versatility  an  examination  of  alternative  turbulator 
shapes  and  geometries  is  possible.  Even  the  evaluation  of  the 
"jersey-barrier”  shape  in  three  dimensions  is  possible  using 
body-fitted  coordinates. 

The  non-linearity  of  the  flow  and  the  large  number  of 
equations  which  are  solved  required  thousands  of  iterations 
for  a  converged  solution.  This  can  take  an  exceptionally  long 
time  even  when  operating  on  a  486-processor  based  computer. 
A  3000  iteration  solution  using  the  configuration  of  this 
report  took  50  hours  to  complete. 

Phoenics  also  requires  the  user  to  accept  certain 

i 

parameters  and  algorithms  without  the  option  to  alter  them. 
Although  numerous  areas  are  accessible,  many  of  the  turbulence 
parameters,  wall  functions  and  other  built-in  subroutines  are 
not.  This  precludes  making  changes  to  fit  the  problem  being 
examined . 

A  full  analysis  of  a  duct  with  any  style  turbulator  should 
include  a  three  dimensional  examination  of  the  flow.  The  500 
X  55  grid  used  here  should  be  considered  to  be  a  minimum  in 
two  dimensions.  The  spanwise  direction  would  add 
significantly  more  points  in  the  grid  and  ultimately  more 
equations  to  be  solved.  This  will  be  excessively  time 
consuming  and  may  or  may  not  provide  the  accuracy  required. 

For  these  reasons  the  recommendation  of  this  report  is  to 
*  use  Phoenics  only  for  the  first  approximations  of  the  solution 
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to  the  duct  problem  using  coarse  grids.  The  ease  with  which 
programs  can  be  entered,  changed  and  run  makes  it  an  excellent 
candidate  for  this  type  of  initial  investigation.  It  can 
locate  areas  of  more  increased  flow  activity  and  aid  in  the 
delineation  of  zones  requiring  clustered  gridding.  Several 
iterations  with  increasingly  finer  and  more  clustered  grids 
will  assist  in  either  supporting  or  refuting  the  initial  grid 
refinements . 

Final  analysis  should  be  carried  out  using  a  program 
possessing  at  least  as  much  power  as  Phoenics  but  with  more 
flexibility  and  a  faster  running  environment.  Two  such 
programs  are  Proteus  3-D  and  Overflow. 
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